o 



A closest vector problem arising in 
radiation therapy planning* 

Celine Engelbeenj Samuel Fiorinif Antje Kiesel^ 
March 12, 2010 

Abstract 



t-h In this paper we consider the following closest vector problem. We are given a set 

of 0-1 vectors, the generators, an integer vector, the target vector, and a nonnegative 

i— i integer C. Among all vectors that can be written as nonnegative integer linear 

combinations of the generators, we seek a vector whose ^-distance to the target 

Q vector does not exceed C, and whose ^-distance to the target vector is minimum. 

^ First, we observe that the problem can be solved in polynomial time provided the 

generators form a totally unimodular matrix. Second, we prove that this problem 
is NP-hard to approximate within an 0(d) additive error, where d denotes the di- 
mension of the ambient vector space. Third, we obtain a polynomial time algorithm 
that either proves that the given instance has no feasible solution, or returns a vector 

CO whose ^-distance to the target vector is within an 0(Vd Ind) additive error of C 

and whose ^i-distance to the target vector is within an 0(dy/d Ind) additive error 
of the optimum. This is achieved by randomly rounding an optimal solution to a 
natural LP relaxation. 

The closest vector problem arises in the elaboration of radiation therapy plans. 
In this context, the target is a nonnegative integer matrix and the generators are 
certain 0-1 matrices whose rows satisfy the consecutive ones property. Here we begin 
by considering the version of the problem in which the set of generators comprises 
?_i all those matrices that have on each nonzero row a number of ones that is at least a 

certain constant. This set of generators encodes the so-called minimum separation 
constraint. We conclude by giving further results on the approximability of the 
problem in the context of radiation therapy. 
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1 Introduction 




Figure 1: A complete treatment unit at Saint-Luc Hospital (Brussels, Belgium). 

Nowadays, radiation therapy is one of the most used methods for cancer treatment. The 
aim is to destroy the cancerous tumor by exposing it to radiation while trying to preserve 
the normal tissues and the healthy organs located in the radiation field. Radiation is 
commonly delivered by a linear accelerator (see Figure [T| whose arm is capable of doing 
a complete circle around the patient in order to allow different directions of radiation. In 
intensity modulated radiation therapy (IMRT) a multileaf collimator (MLC, see Figure [2]) 
is used to modulate the radiation beam. This has the effect that radiation can be delivered 
in a more precise way by forming differently shaped fields and hence improves the quality 
of the treatment. 

After the physician has diagnosed the cancer and has located the tumor as well as 
the organs located in the radiation field, the planning of the radiation therapy sessions is 
determined in three steps. 

In the first step, a number of appropriate beam directions from which radiation will 
be delivered is chosen [TT] . 

Secondly, the intensity function for each direction is determined. This function is 
encoded as an integer matrix in which each entry represents an elementary part of the 
radiation beam (called a bixel). The value of each entry is the intensity of the radiation 
that we want to send through the corresponding bixel. 

Finally, the intensity matrix is segmented since the linear accelerator can only send 
a uniform radiation. This segmentation step mathematically consists in decomposing an 
mxn intensity matrix (or Ruence matrix) A into a nonnegative integer linear combination 
of certain binary matrices S = (sy) that satisfy the consecutive ones property. A vector 
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Figure 2: The multileaf collimator (MLC). 



v G {0, l} d has the consecutive ones property, if v# — 1 and v r = 1 for £ ^ r imply t> j = 1 
for all £ ^ j ^ r. A binary matrix 5* has the consecutive ones property, if each row of 5* 
has the consecutive ones property. Such a binary matrix is called a segment. 
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Figure 3: Leaf positions and irradiation times determining an exact decomposition of a 
4x4 intensity matrix. 

In this paper, we focus on the case where the MLC is used in the so-called step- 
and-shoot mode, in which the patient is only irradiated when the leaves are not moving. 
Actually, segments are generated by the MLC and the segmentation step amounts to 
finding a sequence of MLC positions (see Figure [3]). The generated intensity modulated 
field is just a superposition of homogeneous fields shaped by the MLC. 

Throughout the paper, [k] denotes the set {1,2, ... ,k} for a positive integer k, and 
[£, r] denotes the set {£, £ + 1, . . . , r} for positive integers £ and r with £ ^ r. We also 
allow £ = r + 1 where [£, r] = 0. Thus, an m x n matrix S = (sy) is a segment if and only 
if there are integer intervals [£i, ri\ for i e [to] such that 




1 if je[£i,ri\, 
otherwise. 



A segmentation of the intensity matrix A is a decomposition 

k 

where Uj G Z + and Sj is a segment for j G [fc]. The coefficients are required to be integers 
because in practice the radiation can only be delivered for times that are multiples of a 
given unit time, called a monitor unit. In clinical applications, a lot of constraints may 
arise that reduce the number of deliverable segments. For some technical or dosimetric 
reasons we might look for decompositions where only a subset of all segments is allowed. 
Those subsets might be explicitely given or defined by constraints like the interleaf colli- 
sion constraint (denoted by ICC, also called interleaf motion constraint or interdigitation 
constraint, see [SJ, [UJ, [UJ and [T7j), the interleaf distance constraint (IDC, see [12]). the 
tongue-and-groove constraint (TGC, see [7], [18], [IS], [2D], [21] and [22]), the minimum 
field size constraint (MFC, see [22]), or the minimum separation constraint (MSC, see 

[IS]). 

For some of those constraints, it is still possible to decompose A exactly using the 
set of feasible segments (like for the ICC, the IDC and the TGC), for others an exact 
decomposition might be impossible. In this last case, if S := {Si, . . . , Sk} is our set of 
feasible segments, our aim is to find an approximation B that is decomposable with the 
segments in S, that satisfies 

\\A - BWoo := max - 6y| < C, (1) 

i6[m], j&[n] 

for some given nonnegative integer constant C (possibly, such a matrix B does not exist), 
and minimizes 

\\A-B\\i;= K-M- ( 2 ) 

ie[m],je[n] 

The constraint Q aims at avoiding large bixel-wise differences between target fluence A 
and realized fluence B (that might lead to undesirable hot spots in the treatment), and 
the objective (j2l) measures the total change in fluence with respect to the intensity matrix. 

Later on in this paper, we will focus on the minimum separation constraint, that im- 
poses a minimum leaf opening A G [n] in each open row of the irradiation field. More 
formally, a segment S given by its leaf positions ([£i,ri] , . . . , [£ m ,r m ]) satisfies the mini- 
mum separation constraint if and only if r\ ^ li implies t% — 1% ^ A — 1 for all i G [m\. 
This constraint was first introduced by Kamath, Sahni, Li, Palta and Ranka in [IS], where 
the problem of determining if it is possible to decompose A or not under the minimum 
separation constraint was solved. Here we show that the approximation problem under 
the minimum separation constraint can be solved in polynomial time with a minimum 
cost flow algorithm. 

The approximation problem described above motivates the definition of the following 
Closest Vector Problem (CVP). Recall that the ioo- and li-norms of a vector x G IR d are 
respectively defined by := maxjg^j \xt\ and Hx^ := J2t=i \ x i\- We say that x is 

binary if Xi G {0, 1} for all i G [d]. The CVP is stated as follows: 
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Input: A collection Q = {g 1; g 2 , . . . , g&} of binary vectors in {0, l} d (the generators), a 
vector a in Zl (the target vector), and an upper bound C in Z + U {oo}. 

Goal: Among all vectors b := ^TjLi u jSj with Uj G Z + for j G [k], find one satisfying 
|| a — b || oo ^ C and furthermore minimizing ||a — b||i. If all such vectors b satisfy 
II a — b||oo > C , report that the instance is infeasible. 

Measure: The total change TC := ||a — b||i. 

We remark that the CVP that is the focus of the present paper differs significantly from 
the intensively studied CVP on a lattice that is used in cryptography (see, for instance, 
the recent survey by Micciancio and Regev |24j). 

In order to cope with the NP-hardness of the CVP, we design (polynomial-time, bi- 
criteria) approximation algorithms. For the version of the CVP studied here it is natural 
to consider approximation algorithms with additive approximation guarantees. 

We say that a polynomial-time algorithm is a (A^, A^-approximation algorithm for 
the CVP if it either proves that the given instance has no feasible solution, or returns 
a vector b = X^=i w j§j with Uj G Z + for j G [k] such that ||a — b||oo ^ C + Aoo and 
|| a — b||i ^ OPT + Ai, where OPT is the cost of an optimal solution^ Notice that we 
cannot expect such an approximation algorithm to always either prove that the given 
instance is infeasible or return a feasible solution, because deciding whether an instance 
is feasible or not is NP-complete (this claim holds even when C is a small constant). 

The rest of the paper is organized as follows: Section 2 is devoted to general results 
on the CVP. We start by observing that the particular case where the generators form 
a totally unimodular matrix is solvable in polynomial time. We also provide a direct 
reduction to minimum cost flow when the generators have the consecutive ones property. 
We afterwards show that, when Q is a general set of generators, for all £ > 0, the CVP 
admits no polynomial-time (A^, Ai)-approximation algorithm with A\ ^ (In 2 — e)d, 
unless P = NP. (This in particular implies that the CVP is NP-hard.) We conclude the 
section with an analysis of a natural (Aoo, A ^-approximation algorithm for the problem 
based on randomized rounding [23], with Aoo = 0(\/ d In d) and Ai = 0(dy/d \nd). 

In Section 3, we focus on the particular instances of the CVP arising in IMRT, as 
described above. We first show, using results of Section 2, that the problem can be 
solved in polynomial time when the set of generators encodes the minimum separation 
constraint. We conclude the section with a further hardness of approximation result in 
case A is a 2 x n matrix: for some e > 0, the problem has no polynomial-time (Aoo, Ai)- 
approximation algorithm with Ai ^ en, unless P = NP. (Again, this in particular implies 
that the corresponding restriction of the CVP is NP-hard.) 

In Section 4, we generalize our results to the case where one does not only want to 
minimize the total change, but a combination of the total change and the sum of the 
coefficients Uj for j G [k]. In the IMRT context, this sum represents the total time during 
which the patient is irradiated, called the beam-on time. It is desirable to minimize the 

1 If the instance is infeasible, then we let OPT = oo. 
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beam-on time, for instance, in order to reduce the side effects caused by diffusion of the 
radiation as much as possible. 

Finally, in Section 5, we conclude the paper with some open problems. 



2 The closest vector problem 

In this section we consider the CVP in its most general form. We first consider the 
particular case where the binary matrix formed by the generators is totally unimodular 
and prove that the CVP is polynomial in this case. We afterwards prove that, for all e > 0, 
there exists no polynomial-time (Ac*,, Ai)-approximation algorithm for the general case 
with Ai ^ (In 2 — e) d unless P = NP. We conclude the section by providing a (A^, Ai)- 
approximation algorithm for the CVP, with Aoo = 0(y/d \nd) and Ai = 0(dy/d lad). 

2.1 Polynomial case 

Consider the following natural LP relaxation of the CVP: 

d 

(LP) minJ^ + A) 

k 

s.t. ^ n, • 1, = ai Wie[d] (3) 

3=1 

^ ^ \/i G [d] (4) 

Pi ^ V? G [d] (5) 

ai < C Vi G [d] (6) 

Pi ^ C Wield] (7) 

Uj ^0 Vj G [fc] . (8) 

In this relaxation, the vectors ol and (3 model the deviation between the vector b : = 
X]j=i u jSj an d the target vector a. In the IMRT context, a and f3 model the positive and 
negative differences between realized ffuence and target fluence. Clearly, an IP formulation 
of the CVP can be obtained from (LP) by adding the integrality constraints Uj G Z + for 

3 e >]• 

Let G denote the d x k binary matrix whose columns are gi, g 2 , . . . , gfe. If G is totally 
unimodular, then the same holds for the constraint matrix of (LP). Because a and C are 
integer, any basic feasible solution of (LP) is integer. Thus, solving the CVP amounts to 
solving (LP) when G is totally unimodular. Hence, we obtain the following easy result. 

Theorem 1. The CVP restricted to instances such that the generators form a totally 
unimodular matrix can be solved in polynomial time. 
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2.2 Minimum cost flow problem 

In this section, we assume that the generators satisfy the consecutive ones property. In 
particular, G is totally unimodular. This case is of special interest, because it corresponds 
to the one row case of the segmentation problem in the IMRT context. We show that it 
is not necessary to solve an LP and provide a direct reduction to the minimum cost flow 
problem. 

We begin by appending a row of zeroes to the matrix G and vector a. Similarly, we 
add an extra row to the vectors a and f3. Thus the matrix and the vectors now have d+ 1 
rows. Next, we replace ^ by an equivalent set of equations: We keep the first equation, 
and replace each other equation by the difference between this equation and the previous 
one. Because the resulting constraint matrix is the incidence matrix of a network, we 
conclude that (LP) actually models a minimum cost network flow problem. We give more 
details below. 

We denote the generators by g^ r where [£, r] is the interval of ones of this generator. 
That is, g = ge tr if and only if $ = 1 for i G [£, r] and gi = otherwise. Let I be the set 
of intervals such that Q = {ge,r | [£,r] G X}. We assume that there is no generator with 
an empty interval of ones (that is, I ^ r always holds). Now, let D be the network whose 
set of nodes and set of arcs are respectively defined as: 

V(D) := [d + 1] = {l,2,...,d+l}, and 

A(D) := + | ie[d]}u{(i + l,i) \ie[d)}u{(£,r + l) \ [£,r] el}. 

Let us notice that parallel arcs can appear when the interval of a generator only contains 
one element. In such a case, we keep both arcs: the one representing the generator and 
the other one. 



Figure 4: The network for an instance with d = 6 and k = 9. 

Letting ao := 0, we define the demand of each node j G V(D) as — aj. The arcs 
of type (j, j + 1) and (j + 1, j) have capacity C and cost 1. The other arcs, that is, the 
arcs corresponding to the generators, have infinite capacity and cost 0. An example of 
the network is shown in Figure |4| If we consider a flow p in the network, we have the 
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following correspondence between the flow values and the variables of the LP: 



<p{£,r + l) 
tp(i,i + 1) 



u^ r for all [£, r] G X, 
fa for all iG [of], 
en, for all i £ [d]. 



From the discussion above, we obtain the following result. 

Proposition 2. Let Q and D be as above, let a G and Zei C G Z + U {oo} ; and let 
OPT denote the optimal value of the corresponding CVP instance. Then, OPT equals the 
minimum cost of a flow in D . 

Our network D is similar to the network used in pQ for finding exact unconstrained 
decompositions. There, the arcs of type + 1) and (i + modeling the total change 
are missing and the arcs of type (£,r + 1) are available for all nonempty intervals [£,r]. 

2.3 Hardness 

In this subsection we prove that the CVP is NP-hard to approximate within an additive 
error of at most (In 2 — e)d, for all e > 0. To prove this, we consider the particular case 
where a is the all-one vector. The given set Q is formed of k binary vectors gi, g2, . . . , g*-. 
Because a is binary, the associated coefficients Uj for j G [k] can be assumed to be binary 
as well. 

For our hardness results, we need a special type of satisfiability problem. A 3SAT-6 
formula is a conjunctive normal form (CNF formula) in which every clause contains exactly 
three literals, every literal appears in exactly three clauses and a variable appears at most 
once in each clause. This means that each variable appears three times negated and three 
times unnegated. Such a formula is said to be 5-satisfiable if at most a 5-fraction of its 
clauses are satisfiable. 

As noted by Feige, Lovasz and Tetali [H] , the following result is a consequence of the 
PCP theorem (see Arora, Lund, Motwani, Sudan and Szegedy [3]). 

Theorem 3 ([I!]). There is some < 5 < 1, such that it is NP-hard to distinguish 
between a satisfiable 3SAT-6 formula and one which is 5-satisfiable. 

By combining the above theorem and a reduction due to Feige [13] one gets the following 
result (see Feige, Lovasz and Tetali [H] and also Cardinal, Fiorini and Joret [S]). 

Lemma 4 ((51 EH]). For any given constants c > and C, > 0, there is a polynomial time 
reduction associating to any 3SAT-6 formula $ a corresponding set system <S($) = (V, 5^) 
with the following properties: 

• The sets of 5? all have the size d/t, where d = \V\ and t can be assumed to be 
arbitrarily large. 

• If $ is satisfiable, then V can be covered by t disjoint sets of 5? . 



8 



• If&is 5-satisfiable, then every x sets chosen from 5? cover at most a 1— (l — f) +£ 
fraction of the points, for 1 ^ x ^ ct. 

Theorem 5. For all e > 0, there exists no polynomial-time (A^, Ai)- approximation 
algorithm for the CVP with Ax < (In 2 - e) d « (0.693 - e)d, unless P = NP. 

Proof. We define a reduction from 3SAT-6 to the CVP. We use Lemma [4] to obtain a 
reduction from 3SAT-6 to CVP (by identifying subsets with their characteristic binary 
vectors) with the following properties. For any given constants c > and £ > 0, it is 
possible to set the values of the parameters of the reduction in such a way that: 

• The generators from Q have all the same number | of ones, where t can be assumed 
to be larger than any given constant. 

• If the 3SAT-6 formula <£> is satisfiable, then a can be exactly decomposed as a sum 
of t generators of Q. 

• If the 3SAT-6 formula $ is 5-satisfiable, then the support of any linear combination 
of x generators chosen from Q is of size at most d (1 - (1 - \) x + i) , for 1 < x < ct. 

From what precedes, if $ is satisfiable then the CVP instance is feasible and OPT = 0. 
We claim that if $ is 5-satisfiable, then any approximation b := ^j=i u jSj with Uj G Z + 
for j G [k] has total change TC := | |a — b| |i > d (In 2 — e), provided t is large enough and 
£ is small enough (this is proved below). 

The claim implies the theorem, for the following reason. Assume there exists a 
polynomial-time (A^, Ai)-approximation algorithm with Ai ^ (In 2 — e)d for the CVP 
with some nonnegative integer bound C . Moreover, assume that we are given a 3SAT-6 
formula that is either satisfiable or 5-satisfiable. 

The approximation algorithm either declares the instance given by the reduction to 
be infeasible or provides an approximation b. In the first case, we can conclude that 
$ is not satisfiable, hence 5-satisfiable. In the latter case, we compare the total change 
TC of the solution returned by the algorithm to (In 2 — e) d. If TC ^ (In 2 — e) d then 
the claim implies that $ is satisfiable. If TC > (In 2 — e) d then we can conclude that 
$ is not satisfiable, hence 5-satisfiable, because otherwise the CVP instance would be 
feasible with OPT = and the approximation returned by the algorithm should satisfy 
TC ^ + Ai ^ (In 2 — e) d. In conclusion, we could use the algorithm to decide if $ is 
satisfiable or 5-satisfiable in polynomial time. By Theorem [3j this would imply P = NP, 
contradiction. 

Now, we prove the claim. Notice that we may assume that Uj G {0, 1} for all j G [k]. 
Let x be denote the number of coordinates Uj that are nonzero. We distinguish three 
cases. 

• Case 1: x = 0. 

In this case TC = d > (In 2 — e)d. 
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Case 2 : 1 ^ x ^ ct. 

Let p denote the number of components hi of b that are nonzero. Thus d — p is the 
number of hi equal to 0. The total change of b includes one unit for each component 
of b that is zero and a certain number of units caused by components of b larger 
than one. More precisely, we have: 

TC = d-p+x--p 

t 



> "(f+o-K'-Hy+f) 



x „ / _ 1 



= d \t +2 \ l ~t) - 1 - 26 

= d((l-l3)x + 2/3'-l-2(), 

where f3 :— 1 — |. Note that f3 < 1 and taking t large corresponds to taking /3 close 
to 1. In order to derive the desired lower bound on the total change of b we now 
study the function f(x):—(l — /3)x + 2(5 X . The first derivative of / is 

f(x) = (l-/3) + 21n/3-/r. 

It is easy to verify (since the second derivative of / is always positive) that / is 
convex and attains its minimum at 



ln/3 V. 2111 / 3 
Hence we have, for all x > 0, 

f(x) > /(aw) 

= (l-/3)x min + 2r min 



ln/3 V 2 ln/3 7 ln/3 



By l'Hospital's rule, 



hence we have 



ln/3 V \P-1 



/3->i In /3 



/(a;) >ln2 + l + 2£-e 
for t sufficiently large and £ sufficiently small, which implies 

TC ^ d(ln2 -e) . 
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• Case 3: x > ct. 



Let again p be the number of components hi of b that are nonzero. The first ct 
generators used by the solution have some common nonzero entries. By taking into 
account the penalties caused by components of b larger than one, we have: 




^ d(\n2-e). 

The last inequality holds for t sufficiently large and £ sufficiently small and, for 
instance, c = 2. 

This concludes the proof of the theorem. 

□ 

2.4 Approximation algorithm 

In this subsection we give a polynomial-time (0(y/d\nd), 0(dy/ d\nd )) -approximation 
algorithm for the CVP. This algorithm rounds an optimal solution of the LP relaxation of 
the CVP given in Section 2j] (see page [6]). 



If the LP relaxation (LP) is infeasible, then so is the corresponding CVP instance. 
Now assume that (LP) is feasible and let LP denote the value of an optimal solution of 
(LP). Obviously, we have OPT ^ LP. 

Note that for each basic feasible solution of (LP), there are at most d components of 
u that are nonzero. This is the case, because if we assume that q > d nonzero coefficients 
exist, then only k — q inequalities of type ^ are satisfied with equality. As we have 2d + k 
variables, we need at least 2d + k independent equalities to define a vertex. Thus, there 
must be (2d + k) — (k — q) — d = d + q > 2d independent inequalities of type ([4]), ((5j), ([6]) 
and ([7]) that are satisfied with equality. This is a contradiction, as there can be at most 
2d such inequalities. Thus, for any extremal optimal solution of the linear program, at 
most d of the coefficients Uj are nonzero. 

Algorithm [T] is an application of the randomized rounding technique. This is a widespread 
technique for approximating combinatorial optimization problems, see, e.g., the survey by 
Motwani, Naor and Raghavan [23 J. A basic problem where randomized rounding proves 
useful is the lattice approximation problem: given a binary matrix H of size d x d and 
a rational column vector x G [0, l] d , find a binary vector y G {0, l} d so as to minimize 

Htf(x-y)IL. 

We will use the following result due to Motwani et al. [23] • It is a consequence of the 
Chernoff bound. 

Theorem 6 (|23j). Let (H,x) be an instance of the lattice approximation problem, and 
let y be the binary vector obtained by letting yj = 1 with probability Xj and yj = with 
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Algorithm 1 

Input: a G Z^, C G Z+ U {oo}, and gl , g 2 , . . . , g fc G {0, l} d . 
Output: An approximation b of a. 

If (LP) is infeasible, report that the CVP instance is infeasible. 

Otherwise, compute an extremal optimal solution (ck*,/3*,u*) of (LP). 

for all j G [k] do 

{\u*A with probability u* - [u*\ , 
[u*\ with probability \u*\ - u*. 

end for 

Return b := $^* =1 Ujgj. 



probability 1 — xj, independently, for j G [d]. Then the resulting rounded vector y satisfies 
\\H(pc — y)\\ o0 ^ y/4d\nd, with probability at least 1 — \. 

We resume our discussion of Algorithm [TJ By the discussion above, we know that at 
most d of the components of u* are nonzero. Without loss of generality, we can assume 
that all nonzero components of u* are among its d first components. Then, we let H 
be the d x d matrix formed of the first d columns of G. (W.l.o.g., we may assume that 
d ^ k. If this is not the case we can add generators consisting only of zeros.) Next, we 
let x G [0, l] d be defined via the following equation (where the floor of the u* is computed 
component- wise) : 

<■-*■]-(;). 

Finally, the relationship between the rounded vectors is as follows: 




We obtain the following result. 

Theorem 7. Algorithm^ is a randomized polynomial-time algorithm that either success- 
fully concludes that the given CVP instance is infeasible, or returns a vector b that is a non- 
negative integer linear combination of the generators and satisfies ||a — bH^ ^ C + y/Ad In d 
and || a — b||i ^ OPT + dy/Ad\nd, with probability at least 1 — h. 

Proof. Without loss of generality, assume that (LP) is feasible. Thus Algorithm [l] returns 
an approximation b of a. Let b* = X]j=i M jSi- By Theorem ^ and by the discussion 
above, we have 

llb-bloo = ||G(u-u*) = ||#(x-y) ^ *C VAdlnd, 
with probability at least 1 — \. Now, the result follows from the inequalities 
||a - bHoo ^ ||a - b^loo + ||b* - b||oo < C + ||b* - b||oo 
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and 



la -b||i ^ ||a-b*||i + ||b* - blK < LP + ||b* — blh < OPT + d||b* - b| 



□ 



By a result of Raghavan [26], Algorithm [T] can be derandomized, at the cost of multi- 
plying the additive approximation guarantees v4dlnd and dy/4dhid by a constant. We 
obtain the following result: 



Corollary 8. There exists a polynomial-time (0{y/ d In d ) , 0(dy d In d )) -approximation 
algorithm for the CVP. 

In the case where C = oo, we can slightly improve Theorem [7j as follows. 

Theorem 9. Suppose C = oo. Then, Algorithm^ is a randomized polynomial-time al- 
gorithm that returns a vector b that is a nonnegative integer linear combination of the 

generators and satisfies ||a — b||i ^ OPT + \J 1j y dyfd on average. 

Our proof of Theorem [9] uses the following lemma, which is proved in the appendix. 
Lemma 10. Let q be a positive integer and let Xi,X2, ■ ■ ■ ,X q be q independent random 



variables such that, for all j e [q], P[Xj 



E 



\Xt + X 2 + 



1 - p 3 

+ x q \ 



Pj and P[Xj 



-Pj\ 



Pj . Then 



In 2 



We are now ready to prove the theorem. 
Proof of Theorem [5| We have 



E 



< E 



a 



LP + £ 



LP + E 



E b* 



E 

.i=i 



U i 



Without loss of generality, we may assume that u* = 0, and thus Uj = 0, for j > d. This 
is due to the fact that u* is a basic feasible solution, see the above discussion. 

Now, let Xy := gij (u* — Uj) for all i,j G [d]. For each fixed i 6 [d], Xa, . . . ,X ic i are 
independent random variables satisfying X^ = if g^ = or u* G Z + and otherwise 



X, 



\uj~\ with probability u* — [n*J 
\u*\ with probability \u* \ - w 
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By Lemma 10, we get: 



E 



a — b 



*C LP + E 



£ LP 



+ + • • • + Xid\ 



i=i 



In 2 




d\/d 



<T ;: OPT + dVd. 



LP + J2E\\X a + X t2 + --- + X ld \ 



□ 



A natural question is the following: Is it possible to derandomize Algorithm [T] in order 
to obtain a polynomial-time approximation algorithm for the CVP that provides a total 
change of at most OPT + 0{d\fd) ) provided that C = oo? We leave this question open. 



3 Application to IMRT 

In this section we consider a target matrix A and the set of generators S formed by 
segments (that is, binary matrices whose rows satisfy the consecutive ones property). In 
the first part of this section we consider the case where S is formed by all the segments 
that satisfy the minimum separation constraint. In the last part we consider any set of 
segments S. We show that in this last case the problem is hard to approximate, even if 
the matrix has only two rows. 



3.1 The minimum separation constraint 

In this subsection we consider the CVP under the constraint that the set S of generators 
is formed by all segments that satisfy the minimum separation constraint. Given A e [n], 
this constraint requires that the rows which are not totally closed have a leaf opening 
of at least A. Mathematically, the leaf positions of open rows i G [m] have to satisfy 
Ti — C-i ?3 A — 1. We cannot decompose any matrix A under this constraint. Indeed, the 
following single row matrix cannot be decomposed for A = 3: 

A = ( 1 1 4 1 1 ) . 

The problem of determining if it is possible to decompose a matrix A under this constraint 
was proved to be polynomial by Kamath et al. [T9] . 

Obviously, the minimum separation constraint is a restriction on the leaf openings in 
each single row, but does not affect the combination of leaf openings in different rows. 
Again, more formally, the set of allowed leaf openings in one row i is 

Si = {[£i, Ti] | Ti - £i ^ A - 1 or n = ti - 1}, 
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and does not depend on i. If we denote a segment by the set of its leaf positions 
([£i, ri), . . . , [£ m , r m }) then the set of feasible segments S for the minimum separation con- 
straint is simply <S = «Si X S2 X • • • x S m . Thus, in order to solve the CVP under the 
minimum separation constraint, it is sufficient to focus on single rows. 

Indeed, whenever the set of feasible segments has a structure of the form S = S\ x 
1S2 x • • • x S m , which means that the single row solutions can be combined arbitrarily and 
we always get a feasible segment, solving the single row problem is sufficient. 

From Theorem [TJ we infer our next result. 

Corollary 11. The restriction of the CVP where the vectors are mxn matrices and the 
set of generators is the set of all segments satisfying the minimum separation constraint 
can be solved in polynomial time. 

3.2 Further hardness results 

As we have seen in Subsection 2.1, the CVP with generators satisfying the consecutive 
ones property is polynomial (see Theorem [T]). Moreover, we have proved in Theorem [5] 
that the CVP is hard to approximate within an additive error of (In 2 — e) d for a general 
set of generators. We now prove that, surprisingly, the case where generators contain at 
most two blocks of ones, which corresponds in the IMRT context of having a 2 x n intensity 
matrix A and a set of generators formed by 2 x n segments, is NP-hard to approximate 
within an additive error of en, for some e > 0. 

Theorem 12. There exists some e > such that the CVP, restricted to2xn matrices and 
generators with their ones consecutive on each row, admits no polynomial-time (A^, Ai)- 
approximation algorithm with Ai ^ en, unless P = NP. 

Proof. We prove the theorem again by reducing from the promise problem that was in- 
troduced in Section 12.31 Recall that a 3SAT-6 formula is a CNF formula in which each 
clause contains three literals, and each variable appears non-negated in three clauses, and 
negated in three other clauses. The problem consists in distinguishing between a formula 
that is satisfiable and a formula that is 5-satisfiable. 

Let a 3SAT-6 formula $ in the variables X\, . . . , x s be given. Let c±, . . . , c t denote the 
clauses of $. We say that a variable x« and a clause Cj are incident if Cj involves x« or 
its negation Xj. By double-counting the incidence between variables and clauses, we find 
6s = 3t, that is, t = 2s. 

We build an instance of the restricted CVP as follows: let A = (ay) be the matrix 
with 2 rows and 10s columns defined as 



Thus A has 6s ones, followed by 4s zeros in its first row, and 10s ones in its second row. 
The size of A is 2 x n, where n = 10s = 5t. 




if i = 1 and 6s + 1 < j ^ 10s, 

1 otherwise. 
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To each variable Xi there corresponds an interval I(xi) := [6i — 5, 6i] of 6 ones in the 
first row of A (in this proof, for the sake of simplicity, we identify intervals of the form 
[£,r] and the lxn row vector they represent). We let also I(xi) '■— I(xi). 

For each interval I(xi) = I(xi) = [6i — 5, 6z] we consider two decompositions into three 
sub-intervals that correspond to setting the variable Xi true or false. The decomposition 
corresponding to setting x,- t true is I(xi) = I\{xi) + ^(a^) + h(xi), where Ii(xj) := [6i — 
5, 6i — 5], J 2 (xj) := [6i — 4, 6i — 2] and/ 3 (xj) := [6z— 1, 6«] . The decomposition corresponding 
to setting x^ false is I(xj) = I\{xi) + ^(^i) + ^3(^1), where I\(xi) : = [6i — 5,6i — 4], 
'■— [6i — 3, 6i — 1] and I^Xi) '■— 6i]. An illustration is given in Figure [5j 



-5 -4-3 -2 -1 



h{xi) 



-5 -4-3-2-1 



h{xi) 



h{xi) 



h{xi) 

\ h(Xi) 



Figure 5: The sub-intervals used for the variables. 

Similarly, to each clause Cj there corresponds an interval I(cj) := [5j — 4, 5j] of 5 ones in 
the second row of A. We define ten sub-intervals that can be combined in several ways to 
decompose I{ Cj ). We let h{ Cj ) := [5j-4,5j-4], I 2 ( Cj ) := [5j-2,5j-2], I 3 ( Cj ) := [5j,5j], 
h(cj) ■= [5j-3,5j-3],/ 5 (c,) := [5j-l, 5j-l], I 6 ( Cj ) := [5j-4, 5j-3], J 7 (c,) := [5j-1,5j], 
I 8 (cj) := [5j — 3, 5j — 1], I 9 (cj) := [5j — 4, 5j — 1] and ho(cj) := [5j-3,5j]. An illustration 
is given in Figure |6j 



-4-3-2-1 -4-3-2-1 -4-3-2-1 




Figure 6: The sub-intervals used for the clauses. 

The three first sub-intervals I a {cj) (a £ {1,2,3}) correspond to the three literals of 
the clause Cj. The last seven sub-intervals I a (cj) (a £ {4, ... , 10}) alone are not sufficient 
to decompose I(cj) exactly. In fact, if we prescribe any subset of the first three sub- 
intervals in a decomposition, we can complete the decomposition using some of the last 
seven intervals to an exact decomposition of I(cj) in all cases but one: If none of the three 
first sub-intervals is part of the decomposition, the best we can do is to approximate I(cj) 
by, for instance, Ig(cj), resulting in a total change of 1 for the interval. 
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The CVP instance has k = 20s = lOt allowed segments. The first 3t segments cor- 
respond to pairs (y^Cj) where yi G {xi,Xi} is a literal and Cj is a clause involving y^. 
Consider such a pair (y iy Cj) and assume that yi is the a-th literal of Cj, and Cj is the /3-th 
clause containing y^ (thus a, /3 G {1,2,3}). The segment associated to the pair (yi,Cj) 
is S(yi,Cj) := (Ip(yi),I a { c j))- The last 7i segments are of the form S 1 (cj) := (0, J 7 (c^)) 
where Cj is a clause and 7 G {4, . . . , 10}. We denote the resulting set of segments by 
S = {Si, . . . , Sk}- This concludes the description of the reduction. Note that the reduc- 
tion is clearly polynomial. 

Now suppose we have a truth assignment that satisfies a of the clauses of $. Then we 
can find an approximate decomposition of A with a total change of t — a by summing all 
3s segments of the form S(yi, c,-) where yi = Xi if Xi is set to true, and = Xi if Xi is set 
to false, and a subset of the segments S 7 (cj) for all clauses Cj of the formula. 

Conversely, assume that we have an approximate decomposition B := Ylj=i u j^j °f A 
with a total change of TC := \A — B\\. Because A is binary, we may assume that u is 
also binary. We say that a variable Xi is coherent (w.r.t. B) if the deviation TC(ij) for 
the interval I(x{) is zero. The variable is said to be incoherent otherwise. Similarly, we 
say that a clause Cj is satisfied (again, w.r.t. B) if the deviation TC(cj) for the interval 
I(cj) is zero, and unsatisfied otherwise. 

We can modify the approximation B = Ylj=i u j^j m sucn a wa Y to make all variables 
coherent, without increasing TC, for the following reasons. 

First, while there exist numbers a, i, j and f such that S(xi,Cj) and S(xi,Cj>) are 
both used in the decomposition and the a-th clauses respectively containing 

X{ and Xi, we can remove one of these two segments from the decomposition and change 
the segments of the form Sj(cj) or S^(cji) that are used, in such a way that TC does not 
increase. More precisely, TC(xj) decreases by at least one, and TC(cj) or TC(cy), but 
not both, increases by at most one. Thus we can assume that, for every indices a, i, j 
and j' such that Cj is the a-th clause containing Xi and Cj> is the a-th clause containing 
Xi, at most one of the two segments S(xi,Cj) and S(xi,Cj') is used in the decomposition. 
Similarly, we can assume that, for every indices a, i, j and j' such that cj and cy are the 
a-th clauses respectively containing Xi and Xi, at least one of the two segments S(xi,Cj) 
and S(xi, Cj>) is used in the decomposition. All in all, we have that exactly one of the two 
segments S(xi,Cj) and S(xi,Cj>) is used. 

Second, while some variable Xi remains incoherent, we can replace the segment of 
the form S(yi,Cj), where yi G {xi,Xi}, with a segment of the form S(y~i,Cj/) and change 
some segments of the form S-,(cj) or S 1 (cf) ensuring that TC does not increase (again, 
TC(xj) decreases by at least one and either TC(cj) or TC(cf) increases by at most one). 
Therefore, we can assume that all variables are coherent. 

Now, by interpreting the relevant part of the decomposition of B as a truth assignment, 
we obtain truth values for the variables of $ satisfying at least t — TC of its clauses. 

Letting OPT($) and OPT(v4, S, 00) denote the maximum number of clauses of $ that 
can be simultaneously satisfied and the total change of an optimal solution of the CVP 
with C = 00, we have OPT($) = t - OPT(A, S, 00). 

Using Lemma [4] we know that there exists some 5 G (0, 1) such that it is NP-hard to 
distinguish, among all 3SAT-6 instances $ with t clauses, instances such that OPT($) = t 
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from instances such that OPT($) < 5t. Let e — (1 — S)/5. Using the above reduction 
and the same ideas as in the proof of Theorem |5j the theorem follows. □ 



4 Incorporating the beam-on time into the objective 
function. 

In this section we generalize the results of this paper in the case where we do not only 
want to minimize the total change, but a combination of the total change and the sum of 
the coefficients Uj for j G [k]. More precisely, we replace the original objective function 
||a - b||i by 

k 

fi ■ ||a — b||i + v • Uj, 

i=i 

where [i and v are arbitrary nonnegative importance factors. Throughout this section, we 
study the CVP under this objective function. The resulting problem is denoted CVP-BOT. 

Let us recall that in the IMRT context the generators from Q represent segments that 
can be generated by the MLC The coefficient Uj associated to the segment gj for j G [k] 
gives the total time during which the patient is irradiated with the leaves of the MLC in a 
certain position. Hence, the sum of the coefficients exactly corresponds to the total time 
during which the patient is irradiated (beam-on time). In order to avoid overdosage in 
the healthy tissues due to unavoidable diffusion effects as much as possible, it is desirable 
to take the beam-on time into account in the objective function. 

Here, we observe that the main results of the previous sections still hold with the new 
objective function. 

First, for the hardness results, this is obvious because taking [i = 1 and v = gives 
back the original objective function. 

Second, for showing that CVP-BOT is polynomial when matrix G defined by the 
generators is totally unimodular, we use the following LP relaxation: 



d k 

(LP') minii-J2( a i + + J2 

i=i j=i 



k 



s.t. '^^ u j9ij — oti + & = di Wi G [d] 

a { ^ Vz G [d] 

ft ^ Vz G [d] 

on < C Vz G [d] 

Pi ^ C Vz G [d] 

uj ^0 Vj G [k] . 

Furthermore, if the columns of G satisfy the consecutive ones property, we can still give a 
direct reduction to the minimum cost flow problem. Indeed, it suffices to redefine the cost 



18 



of the arcs of D by letting the cost of arcs of the form + 1) or (j + (for j G [d]) 
be n, and the costs of the other arcs be v. 

Finally, we can also find an (o ( \/ d In d) , O (d \/ d In d) j -approximation algorithm for 

CVP-BOT, by using an extension of the randomized rounding technique due to Srivinasan 
|27j . and its recent derandomization by Doerr and Wahlstrom [9]. 



Consider an instance (H, x) of the lattice approximation problem. Assume that 



Xi g 



Z + . We wish to round x to a binary vector y such that £ J=1 x j = £j=i %' an< ^ ll-^( x — 
y)||oo = O(Vdlnd). Srivinasan [27] obtained a randomized polynomial-time algorithm 
achieving this with high probability. A recent result of Doerr and Wahlstrom [9] implies 
the following theorem. 

Theorem 13 ([9]). Let H e {0, l} dxd and x G Q d n[0,l] d such that J2l=i x j € z +- Then > 



a vector y can be computed in time 0(d 2 ) such that J2j=iVj 



=1 Xj and 



l^(x-y)|| 00 ^(e-l)v / aTnl. 



Let again (ck*,/3*,u*) denote any extremal optimal solution of (LP'). Recall that at 
most d of the k components of u* are nonzero. Without loss of generality, we can assume 
that u* = for j > d. 

yd 



Now, define H and x as previously. Because it might be the case that Y2j=i x j & 



Z 



+ • 



we turn H and x respectively into a (d + 1) x (d + 1) matrix and a (d + 1) x 1 vector by 



letting x d+ i := 



=1 X 3 



By Theorem 
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j=1 and /id+ij 



d+1 



for all i, j G [d + 1]. 



one can find in 0(d 2 ) time a vector y G {0, l} d+1 such that 



(7 ; 



Eji^j and ||iJ(x - y)!!^ ^ (e - l)^(d + 1) ln(d + 1) = O(Vdlnd). We then let 
[u*J + yj for j G [d] and Uj = for j G [A;] \ [<fj. The corresponding approximation of a is 
b := Gu. Notice that the beam-on-time £ J=1 u* will be rounded to |_X)i=i M jJ ^ = -'- 



and to \J2j=i U V\ if Vd+i 



Theorem 



we see that Theorem 



0. Using similar arguments as those used in the proof of 
13 leads to a polynomial-time ( 0(y/d In d), 0(d y/dlnd) 



approximation algorithm for CvP^ 



BOT. 



5 Conclusion 

Here are further questions we leave open for future work: 

• Our results confirm that it is worth to solve the natural LP relaxation of the problem 
(see page [6]). This can be done efficiently when the generators are explicitly given. 
However, in practical applications, the generators are implicitly given. Improving 
our understanding of when the LP relaxation can still be solved in polynomial-time 
is a first interesting open question. 

• Our second question concerns the tightness of the (in)approximability results devel- 
oped here, in particular for the case C = oo. We prove that there is no polynomial- 
time approximation algorithm with an additive approximation guarantee of (In 2 — 
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e) d, unless P = NP. On the other hand, we give a randomized approximation algo- 
rithm with an 0{dy/d) additive approximation guarantee. What is the true approx- 
imability threshold of the CVP (restricted to instances where C = oo)? 

• Our third question is more algorithmic and concerns the application of the CVP 
to IMRT. There exists a simple, direct algorithm for checking whether an intensity 
matrix A can be decomposed exactly under the minimum separation constraint or 
not [T6"t FTP"]. Is there a simple, direct algorithm for approximately decomposing an 
intensity matrix under this constraint as well? 
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Note added during the revision of this manuscript. A recent breakthrough result 
due to Bansal [I] might lead to a randomized (0(yd), 0(dv^))-approximation algorithm 
for the CVP, hence improving both Theorem [7] and Theorem [9j 
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Appendix 



Lemma. Let q be a positive integer and let X\, X%, . . . ,X q be q independent random 
variables such that, for all j G [q], P[Xj = 1 — Pj] = Pj and P[Xj = —pj] = 1 — Pj. Then 



E 



\X 1 + X 2 + --- + X q 



In 2 



y/q- 



Proof. Let X := X\ + X 2 + • • • + X q . By Lemmas A. 1.6 and A. 1.8 from Alon and Spencer 
[2], we have 

E[e xx ] ^ {pe^ x + {l-p)e~ xp X ^ e^ 9 , 
for every A G 1 and p := - Yuj=iPj- Hence, 

e W|] ^ s[e A|x|] ^ E[e x X] + S[e -Axj ^ 2e ^ 

where the first inequality follows from Jensen's inequality (and the fact that the exponen- 
tial is convex) and the second inequality follows from the easy inequality e' x ' ^ e x + e~ x . 

Taking logarithms and letting A = w we get the result. □ 
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